library(tidyverse); theme_set(theme_bw())
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(fs);
library(cowplot); theme_set(theme_cowplot(font_size=6));
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(trend)
# File paths
# roi_file <- "/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_roi/roi_coverage_mapq3.txt"
# wg_path <-"/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_wg"
# af_file <- "/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/Veeramah_PoolSeq_AFs_corrected.filtered_chrXIX.txt"
# out_file <- "/labs/kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/poolseq_depth_figure.pdf"
# File paths
roi_file <- "/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_roi/roi_coverage_mapq3.txt"
wg_path <-"/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/06_samtools_coverage_wg"
af_file <- "/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/Veeramah_PoolSeq_AFs_corrected.filtered_chrXIX.txt"
out_file <- "/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/poolseq_depth_figure.pdf"
out_file_af <-"/Volumes/lab_kingsley/ambenj/myosin_dups/analysis/transplant_pops/gasAcuv5_C4masked/poolseq_AF_figure.pdf"

Depth analysis

# Read roi file
roi <- read_tsv(roi_file) %>% 
  mutate(bam = str_remove_all(bam, "05_mkdup_index/|.PE_ME_merged.sort.Mkdup.realign.BQrecal.realignGA5_C4masked.sort.merged.mkdup.bam|.PE_ME_merged.sort.realign.BQrecal.realignGA5_C4masked.sort.merged.mkdup.bam|.PE_ME_merged.sort.Mkdup.rename.realign.BQrecal.realignGA5_C4masked.sort.merged.mkdup.bam")) %>% 
  rename(samp = bam)
## Rows: 253 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): desc, region, chr, bam
## dbl (8): startpos, endpos, numreads, covbases, coverage, meandepth, meanbase...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
 # separate(bam, into = c("pop", "year"), sep="-", convert = TRUE)
roi
## # A tibble: 253 × 12
##    desc        region chr   startpos endpos numreads covbases coverage meandepth
##    <chr>       <chr>  <chr>    <dbl>  <dbl>    <dbl>    <dbl>    <dbl>     <dbl>
##  1 NLRC5       chrXI… chrX…  2601498 2.63e6    15164    30076     95.3      67.9
##  2 MYH3C1      chrXI… chrX…  2637324 2.65e6     4770    13318     89.2      47.1
##  3 MYH3C2      chrXI… chrX…  2655758 2.67e6     3912    10519     96.8      54.3
##  4 MYH3C3      chrXI… chrX…  2667295 2.68e6     3602    10178     94.9      50.6
##  5 duplicatio… chrXI… chrX…  2665791 2.68e6     5039    16600     94.2      42.0
##  6 SYT19       chrXI… chrX…  2697955 2.73e6    11447    24674     75.5      50.9
##  7 CALB2A      chrXI… chrX…  2745064 2.76e6    13620    14596    100       138. 
##  8 HTRA1A      chrVI… chrVI 14382418 1.44e7    25067    16913     98.7     217. 
##  9 MYH3C3y     chrY:… chrY   8892731 8.90e6     3282    11244    100        43.4
## 10 MYH3C2y     chrY:… chrY   8904789 8.92e6     4912    10952    100        68.2
## # ℹ 243 more rows
## # ℹ 3 more variables: meanbaseq <dbl>, meanmapq <dbl>, samp <chr>
# Read whole genome files
files <- dir_ls(path = wg_path, glob = "*wg_coverage.txt")

# function to add file name to dataframe
read_and_record_filename <- function(filename){read_tsv(filename) %>%
  mutate(filename = path_file(filename))
  }

# gather files into one dataframe
wg <- map_df(files, read_and_record_filename) %>% 
  separate(filename, into=c("samp"), sep="\\.", extra="drop")
wg
## # A tibble: 552 × 10
##    chr   startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq
##    <chr>    <dbl>  <dbl>    <dbl>    <dbl>    <dbl>     <dbl>     <dbl>    <dbl>
##  1 chrI         1 2.96e7 37727036 29260296     98.8      190.      22.4     58.1
##  2 chrII        1 2.37e7 27651994 23455798     99.0      174.      22.5     58.6
##  3 chrI…        1 1.78e7 23018107 17487754     98.5      194.      22.5     57.8
##  4 chrIV        1 3.42e7 41943758 33526396     98.1      183.      22.5     57.9
##  5 chrIX        1 2.08e7 26353055 20426020     98.0      189.      22.5     57.2
##  6 chrM         1 1.57e4    87917    15542     98.7      833.      22.5     47.1
##  7 chrUn        1 1.99e7 21357254 17631961     88.7      162.      22.4     37.3
##  8 chrV         1 1.56e7 21594389 15292087     98.3      208.      22.4     57  
##  9 chrVI        1 1.88e7 22998164 18593175     98.8      183.      22.5     58.5
## 10 chrV…        1 3.08e7 37101356 30240883     98.3      179.      22.4     57.9
## # ℹ 542 more rows
## # ℹ 1 more variable: samp <chr>
wg %>% 
  filter(!chr %in% c("chrM", "chrUn")) %>% 
  ggplot(aes(chr, meandepth)) +
  geom_point() +
  coord_flip() +
  facet_wrap(~samp, nrow = 3)

# get whole genome mean depth (excluding chrUn, chrM, chrXIX, chrY)
wg_cov <- wg %>% 
  filter(!chr %in% c("chrUn", "chrM", "chrXIX", "chrY")) %>% 
  group_by(samp) %>% 
  summarize(weighted_mean_depth = weighted.mean(meandepth, endpos))
wg_cov
## # A tibble: 23 × 2
##    samp    weighted_mean_depth
##    <chr>                 <dbl>
##  1 CH-2009                187.
##  2 CH-2010                342.
##  3 CH-2011                245.
##  4 CH-2012                355.
##  5 CH-2015                263.
##  6 CH-2016                425.
##  7 CH-2017                339.
##  8 LB-1999                107.
##  9 LB-2001                238.
## 10 LB-2002                187.
## # ℹ 13 more rows
# Does it look like even depth ratios? Would expect chrY / chrXIX depth = 0.5   3x for every Y with 50/50 so chrXIX / chrY depth of 3 -> might later years more female skewed...can I correct for this?
sex_ratio <- wg %>% 
  filter(chr %in% c("chrXIX", "chrY")) %>% 
  pivot_wider(id_cols = samp, names_from = chr, values_from = meandepth) %>% 
  mutate(sex_ratio= chrXIX / chrY)
sex_ratio
## # A tibble: 23 × 4
##    samp    chrXIX  chrY sex_ratio
##    <chr>    <dbl> <dbl>     <dbl>
##  1 CH-2009  140.   59.7      2.35
##  2 CH-2010  278.   92.5      3.00
##  3 CH-2011  196.   70.0      2.79
##  4 CH-2012  297.   86.1      3.45
##  5 CH-2015  212.   71.5      2.96
##  6 CH-2016  355.  100.       3.55
##  7 CH-2017  283.   83.0      3.41
##  8 LB-1999   92.4  24.8      3.73
##  9 LB-2001  205.   59.7      3.42
## 10 LB-2002  154.   47.8      3.22
## # ℹ 13 more rows
coverage <- left_join(wg_cov, sex_ratio)
## Joining with `by = join_by(samp)`
coverage
## # A tibble: 23 × 5
##    samp    weighted_mean_depth chrXIX  chrY sex_ratio
##    <chr>                 <dbl>  <dbl> <dbl>     <dbl>
##  1 CH-2009                187.  140.   59.7      2.35
##  2 CH-2010                342.  278.   92.5      3.00
##  3 CH-2011                245.  196.   70.0      2.79
##  4 CH-2012                355.  297.   86.1      3.45
##  5 CH-2015                263.  212.   71.5      2.96
##  6 CH-2016                425.  355.  100.       3.55
##  7 CH-2017                339.  283.   83.0      3.41
##  8 LB-1999                107.   92.4  24.8      3.73
##  9 LB-2001                238.  205.   59.7      3.42
## 10 LB-2002                187.  154.   47.8      3.22
## # ℹ 13 more rows
# Add additional depth info to rois
roi_cov <- left_join(roi, coverage) %>% 
  separate(samp, into = c("pop", "year"), sep="-", convert = TRUE) %>%
  mutate(wg_norm_depth = meandepth / weighted_mean_depth,
          chrXIX_norm_depth = meandepth / chrXIX,
          Y_norm_depth = meandepth / chrY,
         yrs_since_found = case_when(pop == "CH" ~ year - 2009,
                                     pop == "SC" ~ year - 2011,
                                     pop == "RS" ~ year - 2009,
                                     pop == "LB" ~ year - 1985),
         pop = case_when(pop == "RS" ~ "RABS",
                         TRUE ~ pop),
         pop = fct_relevel(pop, c("RABS", "SC", "CH", "LB")),
         desc = as.factor(desc),
         desc = fct_relevel(desc, "NLRC5", "SYT19", "CALB2A", "HTRA1A", "MYH3C1", "MYH3C2", "MYH3C3", "MYH3C1y", "MYH3C2y", "MYH3C3y"))
## Joining with `by = join_by(samp)`
# Is it ok to use all of chrXIX for normalization?
roi_cov %>% 
  filter(pop == "CH" & year == "2017" & str_detect(desc, "MYH")) %>% 
  summarise(sum = sum(numreads),
            sum_covbase = sum(covbases))
## # A tibble: 1 × 2
##     sum sum_covbase
##   <dbl>       <dbl>
## 1 94633       65676
wg %>% 
  filter(samp == "CH-2017" & chr == "chrXIX")
## # A tibble: 1 × 10
##   chr    startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq
##   <chr>     <dbl>  <dbl>    <dbl>    <dbl>    <dbl>     <dbl>     <dbl>    <dbl>
## 1 chrXIX        1 2.06e7 38561406 20182037     98.1      283.      26.2     57.7
## # ℹ 1 more variable: samp <chr>
roi_cov
## # A tibble: 253 × 21
##    desc        region chr   startpos endpos numreads covbases coverage meandepth
##    <fct>       <chr>  <chr>    <dbl>  <dbl>    <dbl>    <dbl>    <dbl>     <dbl>
##  1 NLRC5       chrXI… chrX…  2601498 2.63e6    15164    30076     95.3      67.9
##  2 MYH3C1      chrXI… chrX…  2637324 2.65e6     4770    13318     89.2      47.1
##  3 MYH3C2      chrXI… chrX…  2655758 2.67e6     3912    10519     96.8      54.3
##  4 MYH3C3      chrXI… chrX…  2667295 2.68e6     3602    10178     94.9      50.6
##  5 duplicatio… chrXI… chrX…  2665791 2.68e6     5039    16600     94.2      42.0
##  6 SYT19       chrXI… chrX…  2697955 2.73e6    11447    24674     75.5      50.9
##  7 CALB2A      chrXI… chrX…  2745064 2.76e6    13620    14596    100       138. 
##  8 HTRA1A      chrVI… chrVI 14382418 1.44e7    25067    16913     98.7     217. 
##  9 MYH3C3y     chrY:… chrY   8892731 8.90e6     3282    11244    100        43.4
## 10 MYH3C2y     chrY:… chrY   8904789 8.92e6     4912    10952    100        68.2
## # ℹ 243 more rows
## # ℹ 12 more variables: meanbaseq <dbl>, meanmapq <dbl>, pop <fct>, year <int>,
## #   weighted_mean_depth <dbl>, chrXIX <dbl>, chrY <dbl>, sex_ratio <dbl>,
## #   wg_norm_depth <dbl>, chrXIX_norm_depth <dbl>, Y_norm_depth <dbl>,
## #   yrs_since_found <dbl>

If the sex ratios are 1F:1M: 3:1 chromosome ratio (chrXIX:chrY) If the sex ratios are 2F:1M: 5:1 chromosome ratio If the sex ratios are 1F:2M: 1:2 chromosome ratio

roi_cov %>% 
  filter(desc == "HTRA1A") %>% 
  ggplot(aes(yrs_since_found, sex_ratio, color = pop)) +
  geom_hline(yintercept = 3, linetype = "dashed") +
  geom_line(linewidth = 1.5) +
  geom_point(size = 3) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  xlab("Years since founding") +
  ylab("Sex chromosome ratio (chrXIX / chrY depth)")

roi_cov %>% 
  filter(desc == "HTRA1A") %>% 
  mutate(chrXIX = chrXIX / weighted_mean_depth,
         chrY = chrY / weighted_mean_depth) %>% 
  pivot_longer(cols = chrXIX:chrY, names_to = "chromosome", values_to="depth_1x") %>% 
  ggplot(aes(yrs_since_found, depth_1x, color = pop, group = interaction(chromosome, pop))) +
  geom_line(aes(linetype = chromosome), linewidth = 1.5) +
  geom_point(size = 3) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  ggtitle("Sex chromsome depths") +
  xlab("Years since founding") +
  ylab("Normalized depth (1x genome coverage)")

roi_cov %>% 
  filter(chr!= "chrY") %>% 
  ggplot(aes(yrs_since_found, wg_norm_depth, color = pop)) +
  geom_line(size = 1) +
  geom_point(size = 1) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  xlab("Years since founding") +
  ylab("Normalized depth (1x genome coverage)") +
  facet_wrap(~desc, nrow = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

roi_cov %>% 
  filter(desc %in% c("NLRC5", "STY19","CALB2A", "HTRA1A", "MYH3C1", "MYH3C2", "MYH3C3")) %>% 
  ggplot(aes(yrs_since_found, wg_norm_depth, color = pop)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  xlab("Years since founding") +
  ylab("Normalized depth (1x genome coverage)") +
  facet_wrap(~desc, nrow = 2)

roi_cov %>% 
  filter(chr!= "chrY") %>% 
  ggplot(aes(yrs_since_found, chrXIX_norm_depth, color = pop)) +
  geom_line(size = 1) +
  geom_point(size = 1.5) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  xlab("Years since founding") +
  ylab("Normalized depth (1x chrXIX coverage)") +
#  ylim(0,1.7) +
  facet_wrap(~desc, nrow = 2)

roi_cov %>% 
  filter(desc == "MYH3C3") %>% 
  ggplot(aes(yrs_since_found, chrXIX_norm_depth, color = pop)) +
  geom_line(size = 1) +
  geom_point(size = 1.5) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  xlab("Years since founding") +
  ylab("Normalized depth (1x chrXIX coverage)")

roi_cov %>% 
  filter(chr== "chrY") %>% 
  ggplot(aes(yrs_since_found, Y_norm_depth, color = pop)) +
  geom_line(size = 1) +
  geom_point(size = 1.5) +
#  scale_color_manual(values=c("#9970ab", "#5aae61")) +
  xlab("Years since founding") +
  ylab("Normalized depth (1x chrY coverage)") +
#  ylim(0,1.7) +
  facet_wrap(~desc, nrow = 1)

## Variant analysis

af <- read_tsv(af_file) %>% 
  filter(pos < 2900000, pos > 2500000, `RS-2009` > 0.9) %>% 
  pivot_longer("RS-2009":"LB-2017", names_to = "pop", values_to = "af") %>% 
  separate(pop, sep = "-", into = c("pop", "year"), convert = TRUE) %>% 
  mutate(yrs_since_found = case_when(pop == "CH" ~ year - 2009,
                                     pop == "SC" ~ year - 2011,
                                     pop == "LB" ~ year - 1985,
                                     pop == "RS" ~ year - 2009),
         pop = case_when(pop == "RS" ~ "RABS",
                         TRUE ~ pop),
         pop = fct_relevel(pop, c("RABS", "SC", "CH", "LB")),
         ref_af = 1-af)
## Rows: 200098 Columns: 27
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr  (3): chrom, ref, alt
## dbl (24): pos, RS-2009, CH-2009, CH-2010, CH-2011, CH-2012, CH-2015, CH-2016...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
af
## # A tibble: 18,722 × 9
##    chrom      pos ref   alt   pop    year    af yrs_since_found  ref_af
##    <chr>    <dbl> <chr> <chr> <fct> <int> <dbl>           <dbl>   <dbl>
##  1 chrXIX 2509805 C     T     RABS   2009 0.964               0 0.0360 
##  2 chrXIX 2509805 C     T     CH     2009 0.978               0 0.0220 
##  3 chrXIX 2509805 C     T     CH     2010 0.981               1 0.0190 
##  4 chrXIX 2509805 C     T     CH     2011 0.942               2 0.0580 
##  5 chrXIX 2509805 C     T     CH     2012 0.993               3 0.00700
##  6 chrXIX 2509805 C     T     CH     2015 0.88                6 0.12   
##  7 chrXIX 2509805 C     T     CH     2016 0.907               7 0.093  
##  8 chrXIX 2509805 C     T     CH     2017 0.832               8 0.168  
##  9 chrXIX 2509805 C     T     SC     2011 0.994               0 0.00600
## 10 chrXIX 2509805 C     T     SC     2012 0.965               1 0.0350 
## # ℹ 18,712 more rows
af %>% 
  filter(pos > 2673873, pos > 2674873, pop!="RABS") %>% 
  ggplot(aes(yrs_since_found, ref_af, color = pop, fill=pop)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "smooth", se = TRUE) +
#  geom_point(data = filter(af, pop == "RABS"), size = 2.5, shape=5, alpha=0.1) +
#  geom_smooth(method="lm") 
#  geom_line(data = filter(af, pop != "RABS"),size = 0.1, alpha=0.1) +
    theme_cowplot(8) +
  panel_border(color="black", size=0.75) +
  theme(legend.position = c(0.8, 0.3)) +
  scale_color_manual("Population", values=c("#cc79a7", "#7a669b", "#56b4e9")) +
  scale_fill_manual("Population", values=c("#cc79a7", "#7a669b", "#56b4e9")) +
  xlab("Years since founding") +
  ylab("FW Allele frequency")
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

af %>% 
  filter(pos > 2673873, pos > 2674873, pop!="RABS") %>% 
  ggplot(aes(yrs_since_found, ref_af, color = pop, group = interaction(pos, pop))) +
#  geom_point(data = filter(af, pop == "RABS"), size = 2.5, shape=5, alpha=0.1) +
  geom_line(data = filter(af, pop != "RABS"),size = 0.025, alpha=0.8) +
  stat_summary(
    fun = median,
    geom = 'line',
    aes(group=pop),
    size=3,
    line_type = "dashed") +
    theme_cowplot(8) +
  panel_border(color="black", size=0.75) +
  theme(legend.position = "none") +
  scale_color_manual("Population", values=c("#cc79a7", "#7a669b", "#56b4e9")) +
  scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0, 0))) +
  scale_x_continuous(limits = c(-.2,32.2), expand = expansion(mult = c(0, 0))) +
  xlab("Years since founding") +
  ylab("FW Allele frequency")
## Warning in stat_summary(fun = median, geom = "line", aes(group = pop), size =
## 3, : Ignoring unknown parameters: `line_type`

Final plots

depth_plot <- roi_cov %>% 
  filter(desc == "MYH3C3") %>% 
  ggplot(aes(yrs_since_found, chrXIX_norm_depth, color = pop)) +
  geom_point(data = filter(roi_cov, desc == "MYH3C3", pop == "RABS"), size = 2.5, shape=5, stroke=1.5) +
  geom_line(data = filter(roi_cov, desc == "MYH3C3", pop != "RABS"),size = 1) +
  geom_point(data = filter(roi_cov, desc == "MYH3C3", pop != "RABS"), size = 1.5) +
    theme_cowplot(6) +
  panel_border(color="black", size=0.75) +
  theme(legend.position = "none") +
  scale_color_manual("Population", values=c("#d73027", "#cc79a7", "#7a669b", "#56b4e9")) +
    scale_y_continuous(limits = c(0,2.2), expand = expansion(mult = c(0, 0))) +
  scale_x_continuous(limits = c(-.5,32.5), expand = expansion(mult = c(0, 0))) +
  xlab("Years since founding") +
  ylab("MYH3C3 normalized depth")
depth_plot

# Mann-Kendall test for CH read depth
mk_CH <- mk.test(pull(filter(roi_cov, pop == "CH", desc=="MYH3C3"), chrXIX_norm_depth), alternative = "greater" )
mk_CH
## 
##  Mann-Kendall trend test
## 
## data:  pull(filter(roi_cov, pop == "CH", desc == "MYH3C3"), chrXIX_norm_depth)
## z = 3.0038, n = 7, p-value = 0.001333
## alternative hypothesis: true S is greater than 0
## sample estimates:
##        S     varS      tau 
## 21.00000 44.33333  1.00000
# Mann-Kendall test for SC read depth
mk_SC <- mk.test(pull(filter(roi_cov, pop == "SC", desc=="MYH3C3"), chrXIX_norm_depth), alternative = "greater")
mk_SC
## 
##  Mann-Kendall trend test
## 
## data:  pull(filter(roi_cov, pop == "SC", desc == "MYH3C3"), chrXIX_norm_depth)
## z = 2.403, n = 7, p-value = 0.00813
## alternative hypothesis: true S is greater than 0
## sample estimates:
##          S       varS        tau 
## 17.0000000 44.3333333  0.8095238
mk_SC$p.value
## [1] 0.008130468
# Mann-Kendall test for LB read depth
mk_LB <- mk.test(pull(filter(roi_cov, pop == "LB", desc=="MYH3C3"), chrXIX_norm_depth), alternative = "greater" )
mk_LB
## 
##  Mann-Kendall trend test
## 
## data:  pull(filter(roi_cov, pop == "LB", desc == "MYH3C3"), chrXIX_norm_depth)
## z = 1.6083, n = 8, p-value = 0.05388
## alternative hypothesis: true S is greater than 0
## sample estimates:
##        S     varS      tau 
## 14.00000 65.33333  0.50000
# Variant frequency plot
af_plot <- af %>% 
  ggplot(aes(yrs_since_found, ref_af, color = pop)) +
  geom_point(data = filter(af, pos == 2744769, pop == "RABS"), size = 2.5, shape=5, stroke=1.5) +
  geom_line(data = filter(af, pos == 2744769, pop != "RABS"),size = 1) +
  geom_point(data = filter(af, pos == 2744769, pop != "RABS"), size = 1.5) +
    theme_cowplot(8) +
  panel_border(color="black", size=0.75) +
  theme(legend.position = c(0.8, 0.2)) +
  scale_color_manual("Population", values=c("#d73027", "#cc79a7", "#7a669b", "#56b4e9")) +
    scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0, 0))) +
  scale_x_continuous(limits = c(-.5,32.5), expand = expansion(mult = c(0, 0))) +
  xlab("Years since founding") +
  ylab("FW allele frequency")
af_plot

# Save allele frequency plot
ggsave(
  out_file_af,
  plot = af_plot,
  scale = 1,
  width = 3,
  height = 2.5,
  units = c("in"),
  dpi = 300,
)
final_plot1 <- plot_grid(depth_plot, af_plot, ncol = 1, rel_heights = c(1, 1), labels = c("B","C"), label_size = 12)
final_plot1

final_plot2 <- plot_grid(NULL, depth_plot, nrow = 1, rel_widths = c(1, 1), labels = c("E","F"), label_size = 12)
final_plot2

ggsave(
  out_file,
  plot = final_plot2,
  scale = 1,
  width = 6.5,
  height = 1.75,
  units = c("in"),
  dpi = 300,
)
final_plot_h <- plot_grid(depth_plot, af_plot, nrow = 1, rel_widths = c(1, 1), labels = c("B","C"), label_size = 12)
final_plot_h